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ABSTRACT 

The equation of state of an ideal Fermi gas is expressed in terms of 
Fermi-Dirac integrals. We give formulae for evaluation the Fermi-Dirac integrals 
of orders y 2 , 3 /2, and 5 / 2 and their derivatives in various limits of non- and 
extreme degeneracy and relativity. We provide tables and a Fortran subroutine 
for numerical evaluation of the integrals and derivatives when a limit does not 
apply. The functions can be evaluated to better than 1% accuracy for any 
temperature and density using these methods. 

Subject headings: Equation of State: electrons, Fermi gas 
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1. Introduction 

The electron constituent of neutron star envelopes exists in various stages of degeneracy, 
from a classical gas on the surface to an extremely degenerate Fermi gas in the interior. 
Convective regimes may exist within the envelope ( Urpin 1981). Whether convection does 
occur depends on the magnitude of the temperature gradient compared to the adiabatic 
gradient, which is most conveniently calculated from the adiabatic index 

r 2 = -J— (l-i) 

1 - l/G 

where 

&)/&)] 

An estimate of the role of convection requires accurate values for the temperature T and 
density p derivatives of the Pressure P and energy density E. The envelope constituents 
are electrons, ions, and radiation. The total equation of state, which we treat elsewhere, 
also contains contributions from the Coulomb correction to the ions and takes into account 
the changes in binding energy and electron density due to the variation of ionization level 
with T and p. 

The electrons make a significant contribution to the equation of state (EOS) in the 
nondegenerate regime, and dominate when they are degenerate. We treat the electrons as 
an ideal Fermi gas of arbitrary degeneracy and relativity. We have previously evaluated 
the nonderivative electron EOS by a scheme which interpolates in a two-dimensional table 
of Fermi-Dirac integrals (eq. |[2-7|| below) or uses asymptotic formulae in degenerate, 
nondegenerate, relativistic, and nonrelativistic limits, where applicable. Those limiting 
formulae are described in Bludman & Van Riper (1977). 

The evaluation scheme for the direct function is not sufficiently accurate for the 



r = -— 
PdT 



p dP 
+ P~dp 
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calculation of the EOS derivatives by direct numerical differences. Such a scheme 
becomes ever more hopeless as the degeneracy increases. The fundamental problem is 
obtaining temperature derivatives of quantities which have a vanishingly small temperature 
dependence. This problem also presented itself in our derivation of asymptotic formulae for 
the degenerate T derivatives. Since Urpin (1981) had suggested the possibility of convection 
in the degenerate regime, we were particularly interested in obtaining an accurate adiabatic 
index there. 

Our method of evaluating the derivatives of the EOS relies on the derivatives of the 
Fermi-Dirac integrals, which are found by a scheme exactly analogous to our method for 
the integrals themselves. We prepare a table of the derivatives by accurate numerical 
integration for intermediate degeneracy and relativity and derived formulae (sometimes 
relying on other tabulated functions) in the various limits. The interpolation can be made 
with either a fast second order method or with more accurate polynomial schemes. 

We present those formulae here. We will also publish our tables of the Fermi-Dirac 
functions, their derivatives, modified Bessel functions of the first and second kind, and some 
Fermi integrals on the Astrophysical Journal CDROM. We will also include FORTRAN 
subroutines in which our method is implemented. This paper serves as documentation 
for both tables and subroutines. As such, we will give a brief review of the casting of the 
EOS in terms of the Fermi-Dirac functions and will also show the limiting formulae for the 
functions along with the corresponding formulae for their derivatives. 

The monikers "Fermi-Dirac" and "Fermi" receive no consistent usage in the literature. 
We reserve the term Fermi-Dirac integral (or function) for the bivariate (temperature 
and degeneracy parameter 77) functions defined in equation ( |2-7| ). In the relativistic and 
nonrelativistic limits, these functions reduce to expressions involving integrals (eq. ||4-1|| 
below) which depend on rj alone; we refer to these latter as integrals as the Fermi functions. 
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Numerous studies of the Fermi-Dirac integrals and methods for their rapid evaluation 
exist in the literature (though few of these explicitly treat the derivatives). An excellent 
general reference is chapter 24 of the book of Cox & Giuli (1968), which describes the 
general theory, formulae for many limiting cases, and tabulations of the Fermi-Dirac 
integrals (based on unpublished work by Terry W. Edwards). The literature contains a 
number of schemes for calculating one or more of the integrals. Kippenhahn & Thomas 
(1964; see also Kippenhahn, et al. 1967) give a power series expansion for the evaluation 
of the thermodynamic quantities n, P, and E which is valid for non- and mildly degenerate 
gases. Divine (1965) gives a method, based on third order rational function approximations 
to Fermi integrals of orders Y 2 , 3 / 2 , 2, 5 / 2 , and 3, which gives n, P, and E to a stated 
accuracy of 0.3% for arbitrary degrees of degeneracy and relativity. Guess (1966) considers 
an set of functions (Q n ) equivalent to the Fermi-Dirac functions, and gives series expansions 
for high and low temperatures and degeneracies, accompanied by a tables for the central 
regime where none of those limits apply. Tooper (1969), considering relativistic gases for 
arbitrary degeneracy, derives series expansions in the non- and extremely-degenerate limits 
and discusses methods for numerical integration in the intermediate regime. Beaudet 
& Tassoul (1971) give simple formulae with which n, P, and E can be evaluated to an 
accuracy of several percent for the relativistic and/or degenerate regimes. Bludman & Van 
Riper (1977) give simple formulae, accurate to 0.5% for the semi-degenerate, relativistic 
and nonrelativistic regimes. Nadyozhin (1974) and Blinnikov & Rudzskii (1988) consider 
the limiting cases of extreme relativity and give a number of series expansions. Eggleton, 
et al. (1973) derived fitting formulae for the n, P, and E of an ideal electron gas for a range 
of T and p covering the regime of arbitrary degeneracy and relativity. Evaluation of the 
thermodynamic quantities with 5^ order formulae (their Table 4) agree with values based 
on our main table numerical integrations to 0.04%. An extension of a 4^ n order method 
of Eggleton, et al. given by Pols et al. (1995, Appendix A), while thermodynamically 
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consistent, agrees with our integrations to 0.3%. 

There is a considerable literature dealing with the Fermi integrals. We do not mention 
any of this work here, except to take note of Antia's (1993) rational function approximations 
for Fermi integrals of several half-integral orders with stated maximum relative error of 

io- 12 . 

In the next § we give the formulae for several thermodynamic quantities of an ideal 
Fermi gas in terms of the Fermi-Dirac integrals. Asymptotic formulae in the degenerate limit 
are presented in § 3, including special treatment necessary for the temperature derivatives of 
the pressure and energy density. The treatments in other asymptotic regions are discussed 
in §4. Section 5 covers the numerical details, including the integration and interpolation 
methods and the accuracy and efficiency of the scheme for various interpolation orders. 



2. Thermodynamics of an Ideal Fermi Gas 
We will work throughout in terms of the dimensionless degeneracy and temperature 

where fj, is the chemical potential and m is the mass of the fermion (this theory is also 
applicable to ideal neutron and proton gases). Constants such as the Boltzman constant k-Q 
have their usual meaning throughout this paper. The gas is degenerate (nondegenerate) 
for i] ^> (?7 0). We will refer to relativity regimes based on the value of f3, with 
the gas being relativistic (nonrelativistic) for f3 >> 1 (/3 <C 1). The gas also becomes 
relativistic at high density when the \i mc 2 or, equivalently, when r]j3 >> 1. In practice, 
we take the degenerate (nondegenerate) regime to be rj > 70 (rj < —30) and the relativistic 
(nonrelativistic) regime to be (3 > 10 4 (f3 < 10~ 6 ). 
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The zero of energy for the particles is chosen so that the thermodynamic potential is 

1 + exp (/i — e 



o = _r/,i>r / ^in 

h 6 



(2-2) 



where p is the momentum, g is the statistical weight, and 



e = yj (mc 2 ) 2 + ipc) 2 — mc 2 



(2-3) 



is the kinetic energy With the energy so defined, ]i does not contain the rest mass. (We do 
not consider antiparticles — positrons — in this work; neutron star envelopes do not encounter 
the high T and low p where e + appear in significant numbers.) 

The number density n, pressure P, and energy density (per volume) E of an ideal 
Fermi gas are 



n 



,3^3 



h 3 



(3* F 1/2 ( v ,(3) + (3F 3/2 (r ] ,/3) 



167rV2m 4 c 5 ^ 



3h 3 



F 3/2 (r],(3) + -(3F 5/2 (r],(3) 



and 



E = 



in\/2m 4 c 5 



h 3 



F 3/2 ( V ,[3) + PF 5/2 ( V ,{3) 



where the Fermi-Dirac integral of order k is defined as 



F k (v,P) 



oo X 



o exp(x — 77) + 1 



dx. 



(2-4) 
(2-5) 

(2-6) 



(2-7) 



The energy and pressure derivatives with respect to T and n are 



dP 

On 



t 3 \ dr] 2 drj 



'dF 1/2 | pdF 3/2 



dr] 



dr] 



dE 
dn 



2 n ( 9F 3/2 ndFs/ 2 

= mc j3\ — + f3^^- 
t \ Of] or] 



'dF l/2 + dF 3/ 2 



dr] 



dr] 



(2-8) 



(2-9) 
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and 



dP 
df 



16tt^2 



3h 



wrc" u R± J " ip I r f- 3/2 



f ^ Fs/2 + 2^ 



7 

4' 



' dF 3/2 + 1 p dFs/2 

dr] 2 cfy 



2 Fl/2 + ^ dp 



/2 



+ -/3F 3/2 + /3 : 



3/2 



9/3 



'«yF 1/2 ,0-9^3/2 



+ 0- 



9?7 



(2-10) 



<9£ 



h 



I 2" 



^3/2,7 , 02^5/2 



Mi 



drj 



drj 



2 Fl/2 + P dp 



/2 



+ -/3F 3/2 + /3 : 



3/2 



9/3 



<9^1/2 ^9-^3/2 



9?7 



+ 0- 



(2-11) 



where the derivatives of the Fermi-Dirac integrals are 



drj 



x 



exp (2=2) + exp (2=2 



(2-12) 



and 



9/3 



oo (l + \P 



rx 



o 4 [exp(:r — 77) + 1] 



dx. 



(2-13) 



3. The Degenerate Limit 

3.1. Temperature Derivatives 

The temperature dependence of P and E in the degenerate regime is vanishingly small; 
obtaining accurate temperature derivatives is accordingly problematic. In particular, the 
computer representation of temperature-independent terms in (1) and (2) lacks sufficient 
resolution to ensure cancellations which should occur. We implement these cancellations 
analytically and use the resulting thermal terms in the degenerate temperature derivatives. 
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When i] > 70, the following are used instead of ([2-1 0| ) and ( 2-11 ): 



dP 
df 



n,ED 3/l 3 
where the terms in the numerator are 



k B f3z (Ti — T 2 + T 3 — T 4 



'dF, 



1/2 



To 



'5„ * 0*3/2 7 P 2 dF b/2 



9Fi 



1/2 



drj 

dF 3 / 2 
drj 



+ P 



dp 

0F 3/2 



<9K 



th 

1/2 



0- 



dF 



dp 
th" 

3/2 



drj 



drj 



+ P- 



dF, 



3/2 



dr/ 



dR 



th 



1/2 



dr/ 



+ P- 



dF 



th' 



3/2 



drj 



F$ + l_ (3 dF[ / h 2 
drj 2 c??7 



dF 

2 F ^ 2 + (3 dp 



(3-la) 



1,^5/2 

2^ c>77 



drj 



drj J 



3 T th 
2 -i/2 



Fi h fo + P 



dp 



dp 

(3-lb) 



and 



dE 
df 



n.ED 



3^ 



fen/?* (T 5 - T 6 + T 7 - T 8 



1/2 



drj 



+ P- 



dF, 



3/2 



(3-2a) 



where 
T 5 = 

T 6 = 

T 7 - 



2^/2 + + 2^5/2 + P - W 



-F 



1/2 



+ P^ 



w + ^ F ' i2+p2d -^p- 



dFi/2 + p 9F 3/2 



dF 



th 



1/2 



drj 



drj drj 



P- 



dF. 



th' 



3/2 



drj 



drj 



drj 



dF 



th 



3/2 



drj 

^pth 

2 3/2 



+ P- 



dF t 



th' 



5/2 



drj 
?th 



9F s/2 7 fh 9 dF< 



tlT 

5/2 



and 

r 8 



<9K 



th 



3/2 



<9?/ 



0- 



9F t 



th 



+ ^3/2 + /3 



tlT 



3/2 



dp 



(3-2b) 
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3.2. Degenerate Fermi-Dirac Functions 

The asymptotic degenerate formulae are given in terms of 

y = VA^l (3-3) 

where 

A = 1 + r)(3 (3-4) 

When the degenerate gas is relativistic, ^ > 1 and y ss rjP >> 1. Similarly, in the 



nonrelativistic limit, rjfi <C 1, y ~ \ / 2r/f3 « 1. For most functions, different formulae are 
used depending the value of 7/. 

The following expressions for the Fermi-Dirac functions and the thermal terms are 
only used in the extreme degenerate limit (rj > 70). In that limit, the first few terms in the 
degenerate expansions are sufficient (additional terms in the expansion may be found in 
Cox & Giuli [1968, chapter 24]). The expressions are 

F M * >-A/ 2 + ^-^ (3-5a) 

F}tv,P) = ^rv-t^A (3-5b) 



and 



F 5/2 (V,P) 



F$(V:P) 



1 vr 2 i 3 + 277/3 

7T 2 i 3 + 277/3 
= 7/2 ■ 



6^ v/^T^p' 



1 r vr 2 a 5 + 37y/5 

7T 2 3 5 + 377/5 
= 7/2 ■ 



6V2' Vz + vP' 



(3-6a) 
(3-6b) 

(3-7a) 
(3-7b) 
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where 



and 



\ [yA - ln(y + A)] 



if y > 0.05 



1/2 



b 3 - r y 5 + hv 1 - M + tihy 11 if v < °- 05 



V - i [yA - ln(y + A)] 



/3/2 



if y > 0.05 



3* 

hj^ 5 - ik 7 + mv 9 - lis?/ 11 if y < °- 05 



/5/2 



|yA(l + |y 2 )-|y3_|l n (y + A) if 2/ > 0.1 



±^7 _ JL 9 , .15. 11 

28 y 36 y ~ 704 if 



(3-8) 
(3-9) 

(3-10) 



if y < 0.1 

The small y expansions are used to avoid loss of accuracy due to strong cancellations in 
nonrelativistic limit. 



3.3. Degenerate ^-Derivatives 

The rj derivatives and the corresponding thermal functions are, for all values of the 
relativity parameter y, 



dF 1/2 ( v ,[3) 
drj 

drj 



-^(2 + 77/3)" 



1 



6 r] 2 (2 + 7](3y 



7Y 



6^^1(2 + ^)1 



(3-1 la) 
(3-1 lb) 



drj 
drj 



la i 

(2 + ?7/3) 2 



7T 2 3 + Qr]j3 + 2r] 2 (3 2 
+ ~6 77 2 (2 + 77/3) 2 



TT 2 3 + 6r]P + 2t/ 2 /3 2 



6^2 775(2 + ?7/3) § 



(3-12a) 
(3-12b) 



drj 



V2 



772(2 + 77/3)2 



7T 2 15 + 2077/3 + 6r7 2 /? 2 
¥ 77 2 (2 + 77/3) 2 



1 + 



(3-13a) 



and 
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4. Asymptotic Limits when Not Extremely Degenerate 
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4.1. Arbitrary Degeneracy 

In the relativistic (/? > 1) and nonrelativistic (f3 <C 1) limits the Fermi-Dirac functions 
can be expressed in terms of the simpler Fermi functions G, which only depend on rj. The 
Fermi integrals are given by 

G '("»=i exp(E-„) + l '* C ^ 



with derivatives 



dG k r°° x k 

~2 dx. (4-2) 



drj Jo | exp [ 



2 



cxp 



We require these functions for orders k — 1, 3 /2, 2, 5 / 2 , 3, and 7 / 2 . We have prepared, 
by numerical integration, tables of G and dG/drj, for each of those 7 orders, on a grid of 
integral values of —30 < rj < 70. The evaluation of G and dG/drj in the following formulae 
is accomplished by interpolation. We discuss the tables in more detail below. 



4-1.1. Arbitrary Degeneracy and NonRelativistic 

In the nonrelativistic limit, /3 < 10~ 6 , the Fermi-Dirac functions and their 77-derivatives 
reduce to the Fermi functions and their derivatives: 

Fkijli P) — Gkirf), k = V 2 , 3 / 2 , 5 / 2 (4-3) 

and 

d V ~ ~^T' ^ " /2 ' /2 ' /2 ' (4 " 4) 

The /^-derivatives 

r G fc+1 (77), fc = V 2 , 3 / 2 , 5 / 2 . (4-5) 



dp 4 

involve the next higher order Fermi function. The highest order derivative oc G7/2 only 
appears in the expressions for the /9-derivatives of P and E, multiplied by f3 « 1. Because 
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G 7 /2 is not otherwise required, we carry the G7/2 table in a separate file for the convenience 
of implementations where setting dF 5 / 2 (r), f3)/d/3 = in the nonrelativistic limit is sufficient. 

4-1.2. Arbitrary Degeneracy and Extremely Relativistic 



In the relativistic limit, /3 > 10 4 , 



F k ( V ,(3) = 2(3 
and 



d(3 



%G kH { V l k = y 2 , 3 / 2 , 5 / 2 (4-6) 



dF k (r,, (5) ^^[pdG kH (rj) 



drj V 2 dr] 



k = y 2 , 3 / 2 , 5 / 2 . (4-7) 



4.2. The NonDegenerate Limit — i] < —30 

In the nondegenerate limit, we make use of the well-known (Chandrasekhar 1939) 
relations among the Fermi-Dirac integrals and the modified Bessel functions of the second 
kind Ki and K n (we shall henceforth not explicitly write the adjective "modified"). The 
Bessel functions are defined by 



and 



where 



POO 

KM) = / cosh(t)exp[-0cosh(t)]dt (4-8) 
Jo 

roo 

K n (4>) = / cosh(2t) exp [-0 cosh(£)l dt, (4-9) 
Jo 



. , . exp{x) + exp(-rr) . . 

cosh(x) = — — (4-10) 



is the hyperbolic cosine and 

rh = 



4>=\. (4-11) 



15 



Defining 



£ = exp (r? + </>) , 



the nondegenerate Fermi-Dirac functions are 

F 3/2 (77,/3)~B0i (K n (<f>) - KW) , 
F 5/2 (r], (3) ~ B^ [2KM) + (3/3 - 2) K n (<f>)} , 
the Fermi-Dirac functions are the same as their 77-derivatives, 



dr] 



*F k (r,,0), k= V 2 , 3 / 2 , 5 / 2 , 



and the /^-derivatives are 



9/5 

0*3/2(77,/?) 



<j>{K II -K I )--K I 



and 



.9/5 

dF 5/2 (rj,P) f 
dp 



B& 



2<f>(K I -K H ) + -(5K I + K n ) 



B<j>*< 



-Ku - 20 (2X 7 + + 40 2 - 10) 



(4-12) 

(4- 13a) 

(4-13b) 
(4-13c) 



(4-14) 



(4- 15a) 
(4-15b) 

(4-15c) 



We make use of a combination of expansions and tabulations to evaluate Kj and K n . 



4-2.1. Chebyshev Series Expansions 

Accurate Chebyshev series expansions exist for the Bessel functions. Tooper (1969) 
gives expansions, valid for > 8, for K and Kf. 



e -<t> 00 ' . /if? \ 



(4-16) 
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where the / in the summation means 

OO ' J 00 

£ a k T k = -a + ^2a k T k . (4-17) 

k=0 1 k=l 

(Tooper's expressions [Tables 13 and 14] for the Chebyshev expansions of K and Kj 
contain an extra factor of ir.) The Chebyshev polynomials T k (x) need not be explicitly 
calculated if the series is evaluated recursively. Starting with 



V N+1 = V N+2 = 0, (4-18) 



successive iterates are given by 



and finally the sum 



K = 2 (j - l) K +1 - K + 2 + 4, k = N,-.. 0, (4-19) 

E4^(^-l)=^(^o-^)- (4-20) 

We use Tooper's coefficients from his Tables 13 (j = 0) and 14 (j = 1), for which N = 14, 
to obtain 

M4>) = ~^ (&8 - 62) , KM) = ~^ ~ . (4-21) 



The function K u is given by 



K n {<t>) = K {<t>) + Ikj(<P). (4-22) 





For 4> < 8, Chebyshev expansions from Clenshaw (1962) are applicable: 



00 



/ 



and 



K M ) = - In (Jj £ «& T a (JJ + £ < T 2t \Zj (4.23a) 



M^-.n(^|' a ST,g) + I-|g a ar 2t g) ,, 23b ) 
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Each of the series in (|4-23b|) is evaluated by recursion relations similar to those used 
in the evaluation of (|4-16|) . Starting with Q4-1SQ and using the Chebyshev coefficients a 2k 
tabulated by Clenshaw for N = 17, the recursion 

2 

1 K+i-K+2 + alk, k = N,---0, j = 00, 01, 10, II (4-24) 



is applied until the sum 



N ' 

a 2k T 2k{x) 

k=0 



2^0 j = 00,01,70,71. 



(4-25) 



4-2.2. NonDegenerate Small- f3 Expansions 



In the nonrelativistic limit, (3 — > 0, a power series asymptotic expansion ( Gradshteyn 
& Ryzhik 1980) in (3 is applicable. After expanding Kj and Kjj and collecting terms, 

dF 1/2 (r),P) v^tt , , /_ 3„ 15 „, 105 „„ 105 



F1/2 (v,P) 



dr) 



drj 



^ exp(,) (l + If} - iV + - ^) , (4.26a) 

2 l K IJ \ 8 f 128 f 1024 M 1024 M / ' v ; 

(4-26b) 



„ , ^ dF 3/2 (r},P) 3v^F . ./ 5„ 35 ^ 9 2345 „„ . 



8' 128' 16384' 



„ , m dF 5/2 (ri,f3) 15^ ../ 7 539 o2 
F 5/2 (v, P) = isr^ - — ?~ ex P(^) ( 1 + oP - TT^/ 3 



8' 4090' 



(4-26c) 



and 



16 PV /; V 8 M 128 M 32 M 



9/3 32 
dF s/2 (n,P) 105^ 



8 2048 
77 



6XP( ^ ( X - 256^) ' 



(4-27a) 
(4-27b) 

(427c) 



dp 34 

The higher order terms in the expansion ( 4-27c| ) have cancelled, but since this expression 
is only used for ft < 10~ 3 , the remaining terms give sufficient accuracy. 
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4-2.3. Extremely Relativistic and NonDegenerate 

In the relativistic limit, which we employ for (3 > 10 4 , Kj(x) ~ Kjj(x) ~ 2/x 2 , 

and ( |4-13a|) through ( f4-15c| ) become 



and 



F k {„, 8) = 9Ft( ^' l3) ~ (* + 1) ! exp + ±) ^ k — V 2 , 3 / 2 , 5 / 2 (4-28) 



dFti "'' 3) F*,,ft, k = >/„ '/„ '/,. (4-29) 



a/3 v 2 / 5 



4-2.4- Bessel Function Tabulation and Evaluation Methods 

We have prepared, by numerical integration, tabulations of the Bessel functions on a 
grid of —4 < log lo < 1.6 with a spacing Alog lo = 0.1. 



For <p > 1, evaluations of derivatives by ( [4- 15a ), ( 4-15b ), and especially ( 4-15c ) involve 



cancellation between the Kj and Kjj terms and require more precision than exists in our 
table. In addition, Kj and Kjj decrease ever more rapidly with increasing this sharp 
falloff causes loss of interpolation accuracy for > 1 (j3 < 1). Greater accuracy is obtained 
by use of high (6 ) order interpolation in Kj and Kjj and the second expression in 
equations |4-15a] , |4-15b| , and |4-15c| , which result in better cancellation, but the tables should 



only be used above (3 = (3 t . Between (3 t and the smaller /3 C , the Chebyshev expansions 
should be used. For (3 < (3 C , the small-/? expansions are more accurate than the Chebyshev 
series, for which the accuracy suffers as (3 decreases. 

th 

For 6 order interpolation, the tables, rather than the Chebyshev series are used 
for (3 > (3 t = 10 03 . Across this boundary, the Fermi-Dirac functions and derivatives 
match to 1 part in 5 x 10 5 or less. For less accurate, lower order interpolations, larger 
values of j3t may be appropriate (unless lower order interpolation is used for the sake of 
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reducing computer time). The computer time taken to evaluate a set of functions with the 

th 

Chebyshev expansions is 1.09 times greater than evaluating the same set with 6 order 
table interpolation. 

When the /3-derivatives are required, the small-/? expansions should be used for 
(3 < (3 C = 10~ 2 - 8 ; for ^l/lxh^l ; the loss of accuracy with increasing or decreasing (3 is 
steep away from the respective side of the switching point. When the derivatives are not 
required, log/3 c = —2.0 is recommended. Across the (3 C = 10~ 2 ' 8 boundary, the functions, 
the ^-derivatives, and 9Fl ^'^ ma tch to a relative difference of 10~ 6 or better, 8F 3/jifa/3) to 
10~ 5 , and 9F5/ Qp'^ to 0.003. The evaluation of the same set of functions with the Chebyshev 
expansion requires 3.3 times as much computer time as with the small-/? expansions. 

An alternative method of calculating the nondegenerate quantities, which we do not 
employ, makes use the values of the Fermi-Dirac functions and derivatives in the central 
table. The r\ dependence is given by equation ^4-12|, so that 



X{r), (3) = exp^ + 30) X(-30, /?), 
where X is any Fk or derivative. 



5. The Tables: Creation, Interpolation, and Accuracy 

5.1. Numerical Integration 

All integrals considered here were numerically evaluated using a 7-point adaptive 
Newton-Cotes quadrature rule (implemented in the routine QNC79 from the SLATEC 
software library). For the Fermi-Dirac and Fermi integrals, r\ + 100 (rather than oo) was 
used for the upper limit of the numerical integration; the lower limit was max(0,?7 — 100), 
whereas for the Bessel functions the integrations ran from to 6 — ln(</>). The integration 
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routine evaluates each integral until a desired accuracy is achieved. The accuracy is 
expressed as a tolerance S, where the result of the numerical integration / does not deviate 
more than SI from the true answer. We used £ = 10~ 12 . A comparison of integrations 
made on 2 platforms — a Silicon Graphics Indigo2 with an R4400 cpu chip running IRIX 5.3 
and an Apollo DN4000 running DOMAIN/OS 10.3 — agreed to within a relative difference 
of 1.2 x 10~ 7 . This comparison suggests the accuracy of our integrations is not better than 
1 part in 10~ 7 . The values in the table are the Silicon Graphics integrations. 



5.2. The Tables 

Three tables of integrals accompany the electronic version of this paper. The 
file f dints. tab contains the Fermi-Dirac integrals and the Fermi integrals. The file 
df dints. tab contains the derivatives of the Fermi-Dirac and of the Fermi integrals. The 
Bessel functions are found in bessel.tab. The file fermi7h.tab contains the Fermi 
integrals of order 7 / 2 . These tables are read by the FORTRAN subroutine calcdfi, 
which is also provided along with the tables; the exact format of the tables can be found 
by examining the relevant READ statements in the subroutine. The functions and the 
derivatives are maintained as separate files to facilitate implementations where only the 
function values are required; subroutine calcf i, also supplied, is a variant of calcdfi with 
the derivatives stripped out. 

Subroutine calcdfi reads the data from the binary (or unformatted) files f dints. unf 
(which contains the Bessel data) and df dints .f is they both exist; if they do not exist, the 
formatted files f dints .tab, • • • are read, and the logarithm of the data are written to the 
binary files. The binary files are preferred because the formatted files are 5 times as large 
and take 35 times longer to load. All data is held in memory as (natural) logarithms. 
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The main table of the Fermi-Dirac functions and their integrals for —30 < rj < 70 and 
—6 < log/3 < 4, the same ranges as used in the tabulations in Appendix A. 2 of Cox & 
Giuli (1968). The table uses a finer grid for small values of \rj\ and \/3\. The table contains 
21 p values corresponding to the integral and half-integral values of —6 < log(/3) < 4. The 
47 points in the rj grid are concentrated towards rj — 0. The spacing is 0.1 for — 1 < rj < 1, 
1 for —5 < i] < — 1 and 1 < rj < 5, and 5 for —30 < i] < —5 and 5 < i] < 70 . The 77 grid for 
the Fermi integrals and derivatives consists of 101 values of integral —30 < r] < 70. The 
grid for the Bessel functions contains 59 values with a spacing A log 10 = 0.1. 



All interpolation is made with the logarithms of the integrals and derivatives as 
functions of log (3 (or log 0) and r] (except as noted in the discussion of second order 
interpolation in the 2-dimensional table). Use of the direct values, rather than the 
logarithms, is much less accurate for all interpolation orders. 



Second order interpolation for the 1-dimensional functions is a simple linear 
interpolation between the values on bracketing grid points. For the central 2-dimensional 
tables, the interpolation is logarithmic in the rj direction 



5.3. Interpolation 



5.3.1. Second Order Interpolation 



f(P) = exp{(l -d)ln[f(r] a ,P)} + dln[f(r] b ,P)}} 



(5-1) 



where d is the interpolation coefficient and i]^ a ^ are points in the table. The | 
used for the interpolation in the f3 direction, 



f power is 



/= \(l-d)f(p a )i + df^} 1 



(5-2) 
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where d and /3( a ,6) have meanings similar to the above. 

5.3.2. Higher Order Polynomial Interpolation 

Our polynomial interpolation is based on the subroutines POLINT and P0LIN2 from 
the Numerical Recipes book ( Press, et al. 1986). We examined the accuracy and speed of 
the interpolation for several polynomial orders. We found 6 th order to be most accurate, 
but with a substantial penalty in execution speed. Orders 2 through 6 are available in the 
subroutine we provide. 

The tables below give the accuracy and execution time for various orders and functions. 
The accuracy is given in terms of the relative error max(////, f f /f), where / is truth, as 
given by a numerical integration, and // is the value from interpolation. Execution times 
are normalized to 1 for 2 nd order. The timings were made with uniform samplings over the 
respective tables. 

5.3.3. The Two Dimensional Tables 

Table [I] gives the maximum relative error for the Fermi-Dirac integrals, along with the 
execution times, and Table |] gives the maximum relative errors for the derivatives. The 
accuracy is based on comparisons with numerical integrations made at rj — 0.6r/j + 0.4r/j +1 
and log P = 0.6 log (3j + 0.4 log (3j + i for each (i, j) cell in the table, a similar set of integrations 
with 0.4 «-> 0.6 for every other (i,j) cell, and a sampling of cells with a set of 10 integrations 
crossing the cell in some direction. 

The order in which the 1-dimensional interpolations were performed (ie. rj direction 
first or (3 direction first) made no difference in the accuracy. 
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EDITOR: PLACE TABLE | HERE. 
EDITOR: PLACE TABLE | HERE. 

5.3.4- Fermi Integrals 

Table |3| gives the maximum error in the Fermi integrals and their derivatives for a 
range of ±0.5 in rj centered on the value listed. There is variation of up to 100% among the 
relative errors for the individual functions and derivatives. The maximum relative errors 
for the 6 functions and the 6 derivatives are the same to within 1%. The relative execution 
times for the interpolations are given in Table |j. 

EDITOR: PLACE TABLE | HERE. 
EDITOR: PLACE TABLE | HERE. 

5.3.5. Bessel Function Interpolation 

Table | lists the maximum relative errors in the Bessel functions for several ranges 
of 0. For each range in each order, the maximum relative errors for Kj and Kjj were 
comparable. The larger of the two is listed in the table. The error becomes increasingly 
worse as <fi increases and the Bessel functions are dominated by the factor exp —(/). 

EDITOR: PLACE TABLE | HERE. 



5. 3. 6. Accuracy at Treatment Boundaries 



Different methods are used to evaluate the Fermi-Dirac integrals and their derivatives 
for different ranges of rj and /3. The closeness of the evaluations on opposite sides of a 
treatment boundary is given in Table |6|. The quantity tabulated is the largest value of 
max(/ r //i, fi/ fr) along a boundary For an r\ boundary /( r n = j{r\ ± 0.000001), and for a 
(3 boundary = /(log/? ± 0.000001). The results are most excellent, especially for 6 th 
order interpolation. 

The /3-derivatives lose all accuracy near log/3 = —1.6. Accordingly, we also list the 
closeness along 77 = — 30 excluding a range near log/3 = —1.6. 

EDITOR: PLACE TABLE | HERE. 

We are grateful to Onno Pols for helping us correctly implement the Eggleton, Faulkner, 
& Flannery formulae, and to Jim Lattimer for discussions on various approximations. 
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Table 1. Main Table Interpolation. 



Maximum Relative Error 



N time F 1/2 ( V ,P) F 3/2 (rj,P) F 5/2 ( V ,P) 



2 


1.0 


1.082 


1.118 


1.143 


3 


1.8 


1.572 


1.646 


1.628 


4 


2.9 


1.023 


1.023 


1.021 


5 


4.6 


1.027 


1.016 


1.009 


6 


7.1 


1.007 


1.003 


1.004 


7 


10. 


1.005 


1.010 


1.012 


8 


14. 


1.013 


1.007 


1.005 


9 


19. 


1.096 


1.050 


1.026 
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Table 2. Main Table Derivative Interpolation. 



Maximum Relative Error 



dF 1/2 (r,,0) dF 3/2 (r,,P) dF 5/2 (r,,0) dF 1/2 (r,,0) dF 3/2 (r,,P) dF s/2 (r,,P) 
dr, d v d v dp dp dp 



2 


1.034 


1.082 


1.118 


1.132 


1.159 


1.177 


3 


1.391 


1.572 


1.646 


1.573 


1.646 


1.628 


4 


1.019 


1.023 


1.024 


1.024 


1.023 


1.018 


5 


1.032 


1.027 


1.016 


1.022 


1.010 


1.010 


6 


1.011 


1.007 


1.003 


1.005 


1.004 


1.005 


7 


1.019 


1.006 


1.010 


1.007 


1.012 


1.013 


8 


1.026 


1.013 


1.007 


1.010 


1.006 


1.005 


9 


1.187 


1.096 


1.050 


1.069 


1.036 


1.019 
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Table 3. Fermi Function Interpolation Relative Errors. 



Interpolation Order 



V 




2 




3 


4 




5 




6 


-29.5 


1 


.0000006 


1 


.0000006 


1.0000006 


1 


.0000006 


1 


.0000006 


-19.5 


1 


.0000006 


1 


.0000006 


1.0000006 


1 


.0000006 


1 


.0000006 


-14.5 


1 


.0000005 


1 


.0000023 


1.0000007 


1 


.0000008 


1 


.0000006 


-9.5 


1 


.0000068 


1 


.0000182 


1.0000011 


1 


.0000006 


1 


.0000005 


-4.5 


1 


.0009817 


1 


.0025651 


1.0001572 


1 


.0000744 


1 


.0000216 


-2.5 


1 


.0063697 


1 


.0168104 


1.0008456 


1 


.0003333 


1 


.0000158 


-1.5 


1 


.0137371 


1 


.0364101 


1.0011945 


1 


.0002281 


1 


.0002464 


-0.5 


1 


.0225276 


1 


.0566500 


1.0005099 


1 


.0009853 


1 


.0002920 


0.5 


1 


.0247931 


1 


.0412589 


1.0026541 


1 


.0018509 


1 


.0004923 


1.5 


1 


.0187399 


1 


.0212393 


1.0022867 


1 


.0005230 


1 


.0003849 


2.5 


1 


.0132760 


1 


.0471996 


1.0007911 


1 


.0016337 


1 


.0002560 


4.5 


1 


.0091014 


1 


.0247539 


1.0006357 


1 


.0002153 


1 


.0000434 


9.5 


1 


.0041356 


1 


.0051802 


1.0000569 


1 


.0000178 


1 


.0000018 


14.5 


1 


.0020823 


1 


.0019885 


1.0000169 


1 


.0000033 


1 


.0000004 


19.5 


1 


.0012197 


1 


.0008866 


1.0000060 


1 


.0000016 


1 


.0000005 


24.5 


1 


.0007939 


1 


.0004615 


1.0000026 


1 


.0000006 


1 


.0000005 


29.5 


1 


.0005556 


1 


.0002668 


1.0000014 


1 


.0000004 


1 


.0000004 


34.5 


1 


.0004097 


1 


.0001662 


1.0000010 


1 


.0000005 


1 


.0000005 


39.5 


1 


.0003149 


1 


.0001134 


1.0000009 


1 


.0000006 


1 


.0000005 


44.5 


1 


.0002490 


1 


.0000782 


1.0000005 


1 


.0000005 


1 


.0000005 


49.5 


1 


.0002017 


1 


.0000569 


1.0000004 


1 


.0000007 


1 


.0000005 


59.5 


1 


.0001400 


1 


.0000308 


1.0000005 


1 


.0000006 


1 


.0000005 


69.5 


1 


.0001032 


1 


.0000216 


1.0000021 


1 


.0000034 


1 


.0000067 
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Table 4. Interpolation Times. 



Relative Time 



N G, dG/dr] K h K n 



2 


1.0 


1.0 


3 


7.6 


4.9 


4 


11. 


7.1 


5 


16. 


9.8 


6 


22. 


13. 
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Table 5. Bessel Function Interpolation Relative Errors. 



Interpolation Order 
log 4> a log b 2 3 4 5 6 



-4.0 -2.0 1.000066 1.000139 1.000042 1.000047 1.000042 

-2.0 0.0 1.005176 1.008385 1.000117 1.000050 1.000041 

0.0 0.5 1.018382 1.028126 1.000306 1.000092 1.000059 

0.5 1.0 1.060710 1.090631 1.000918 1.000201 1.000088 

1.0 1.5 1.205830 1.311856 1.002875 1.000673 1.000389 
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Table 6. Difference at Treatment Boundaries. 



Interpolation Order 



Y J\ ) LI 11LI CI J- V 


TTi i n pf i on 

J. Lill^ulwll 


2 


3 


4 


5 




t] = 70 


F 


1.0007 


1.0026 


1.0022 


1.0582 


1.0309 


r] = 70 


dF/dr] 


1.0008 


1.0024 


1.0030 


1.0632 


1.0308 


7] = 70 


dF/d(3 


1.0011 


1.0023 


1.0038 


1.0611 


1.0279 


7] = -30 


F 


1.0027 a 


1.0022 


1.0029 


1.0473 


1.0303 


77 = -30 


dF/d-q 


1.0027 a 


1.0022 


1.0029 


1.0473 


1.0303 


7] = -30 


dF/df3 


1.002P 


1.0020 


1.0024 


1.0445 


1.0279 


log/5 = 4 


F 


1.0072 


1.0221 


1.0231 


1.5966 


1.1481 


log/5 = 4 


OF/Ot] 


1.0081 


1.0301 


1.0231 


1.5967 


1.1270 


log/5 = 4 


dF/dp 


1.0064 


1.0207 


1.0225 


1.5971 


1.1481 


log (3 = — 6 


F 


1.0071 


1.0262 


1.0227 


1.5976 


1.1386 


log/3 = -6 


OF/Ot] 


1.0101 


1.0312 


1.0227 


1.5752 


1.1132 


log/5 = -6 


dF/d(3 


1.0042 


1.0134 


1.0213 


1.5864 


1.1258 



a The difference is 1.0007 when the range 3.7 < log/5 < 4.0 is 
excluded 
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